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ABSTRACT 


The  in-plane  edge  impact  of  composite  plates,  with  or  without  a  pro¬ 
tection  strip,  is  investigated  in  this  work.  A  computational  analysis 
based  on  the  Fast  Fourier  Transform  technique  is  presented.  The  particular 
application  of  the  present  method  is  in  the  understanding  of  the  foreign 
object  damage  problem  of  composite  fan  blades.  However,  the  method  is 
completely  general,  and  may  be  applied  to  the  study  of  other  stress  wave 
propagation  problems  in  a  half  space.  Results  indicate  that  for  the  pro¬ 
tective  strip  to  be  effective  in  reducing  impact  stresses  in  the  composite 
the  thickness  must  be  equal  or  greater  than  the  impact  contact  dimension. 
Also  large  interface  shear  stresses  at  the  strip  -  composite  boundary  can 
be  induced  under  impact. 
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INTRODUCTION 

This  report  is  part  of  a  continuing  effort  by  NASA  to  understand  the  basic 
mechanics  of  foreign  object  impact  of  composite  materials  of  particular  interest 
are  damage  resistant  designs  of  jet  engine  fan  blades  under  hail  or  bird 

impact.  In  previous  reports  the  central  or  normal  impact  response  of  composite 

plates  was  examined  [1],  In  this  report  the  mechanics  of  edge  impact  of 

composite  plates  are  examined.  This  is  schematically  illustrated  in  Figure  1. 
The  basic  approach  to  the  study  of  impact  of  composite  plates  in  this 

program  has  been  to  examine  the  stress  waves  generated  by  the  impact  forces. 

For  central  impact  of  plates  it  has  been  shown  that  in  addition  to  wave  pro¬ 
pagation  across  the  plate  thickness,  bending  and  extensional  waves  propagate 
away  from  the  impact  site.  The  stresses  associated  with  these  waves  have  been 
studied  without  considering  the  effect  of  boundaries  such  as  the  free  or  clamp- 

led  edges  of  a  fan  blade.  This  simplification  has  been  made  on  the  premise  that 

-4 

for  short  impact  times  e.g.  less  than  10  sec.  few  edge  reflections  have 
taken  place,  and  that  the  highest  stresses  occur  at  the  impact  site.  However, 
for  edge  impact,  the  boundary  conditions  greatly  affect  the  nature  of  the  wave 
mechanics . 

Edge  waves  in  solids  have  been  studied  extensively  in  seismology.  The 
principal  phenomenon  is  the  entrapment  of  wave  energy  in  a  layer  near  the  sur¬ 
face.  This  surface  wave  is  known  as  a  Rayleigh  wave  and  travels  at  a  velocity 
below  the  shear  velocity  for  isotropic  solids.  For  plates,  however,  two  types 
of  edge  waves  can  occur  as  shown  in  Figure  2a.  For  impact  transverse  to  the 
plate,  flexural  edge  waves  can  occur.  For  in-plane,  Rayleigh,  type 
edge  waves  are  generated.  In  this  report  only  in-plane  edge  waves  will  be 
discussed. 
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Wave  type  solutions  to  the  equations  of  elastodynamics  of  an  orthotropic 

plate  which  exponentially  decay  away  from  the  edge  (X  direction)  and  pro- 

3 

pagate  along  the  edge  (X  direction)  can  be  found  provided  the  edge  wave 

1 

velocity  satisfies  the  equation  (Reference  9) . 


1/2  2 

[C  -  C  (C  -  pv2)]  =  0 
13  33  11 


where  P  is  the  massdensity  of  the  plate  and  are  the  effective  plate 
elastic  constants  (denoted  by  in  Reference  1) .  It  can  be  shown  that 

ij 

one  real  root  lies  in  the.  interval 

0  <  pv2  <  C 

55 

Thus,  the  edge  or  Rayleigh  wave  speed  is  less  than  the  shear  speed  in  this 
1/2 

direction  [C  /p] 

55 

Changing  the  layup  angle  will  affect  the  elastic  constants  C  and  hence, 

ij 

change  the  value  of  the  edge  wave  velocity. .  The  results  of  this  calculation 
are  shown  in  Figure  2b  where  the  C  are  obtained  from  Reference  7  .  The 

ij 

edge  wave  speed  seems  to  obtain  a  maximum  between  _+  15  and  +_  30°  layup  angle  which 


A 


is  below  the  extentional  wave  speeds  labelled  "dilational"  and  "shear" 
in  Figure  2b  and  which  is  greater  than  the  bending  wave  velocity. 
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In  order  to  prevent  damage  to  composite  fan  blades  under  foreign  object 
impact  leading  edge  protection  has  been  used.  (See  e.g.  Ref.  (8)).  This  usually 
consists  of  a  strip  of  metal  attached  to  the  leading  edge  of  the  fan  blade.  To 
model  the  effects  of  this  impact  protection  strip,  the  in-plane  edge  impact 
of  an  anisotropic  plate,  with  a  beam- strip  attached  to  the  impact  edge,  has 
been  studied  (see  Figure  1).  It  will  be  shown  later  that  the  strip  will  decrease  the 
tensile  stress  along  the  edge  while  producing  shear  stress  between  the  strip 
and  the  plate  edge. 
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I.  BASIC  EQUATIONS  FOR  EDGE  LOADING  ON  ANISOTROPIC  PLATE 

Let  the  plate  be  the  half-space  x  j  >  0.  The  equations  of  motion 
are  given  by  [1]  as: 


C  u  +  C  u  +  (C  +C)u  =  p  U] 

11  1,11  55  1,33  13  55  3,13  ,tt 


(1.1) 


C  u  +C  u  +  (C  +C)u  =  p  u  .  . 

33  3,33  55  3,11  13  55  1,13  3,tt 


here  we  have  employed  C  ,  C  ,  C  ...  to  denote  C  ,  C  ,  C  ...  of 

11  13  33  11  13  33 

[1] .  The  in-plane  motion  is  assumed  to  be  independent  of  the  bending  defor¬ 
mation.  With  the  loading  condition  shown  in  Figure  1,  the  boundary  conditions 
are  (without  protection  strip) : 


t  (x  ,0,t)  =  C  (u  +  u 

131  551,3  3 


=  0 

x  =0 

3 


(1-2) 


t  (x  ,0,t) 

3  3  1 


The  following  nondimens ional  parameters  are  used; 
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c*  =  C  /C  ,  C*  =  C  /C  ,  C*  =  C  /C  ,  C*  =  C  /C  ,  P*  =  p/C  / 

11  11  66  13  13  66  33  33  66  55  55  66  0  66 

(1.3) 

u*  =  u  /l,  u*  =  u  n,  x*  =  x  /l,  x*  =  X  /£,  t*  =  t  v^T'/pH  (1.4) 


£  is  a  length  parameter.  C  is  a  typical  elastic  constant.  In  what 

6  6 

follows  the  equations  will  be  assumed  to  be  nondimens ionali zed  using  (1.3) 
and  (1.4). 

To  obtain  the  solution,  we  employ  transform  methods.  Define. 


F(f)  = 


e  lkiXi  f(xi)  dxi 


the  Fourier  transform  of  f(x)  where  we  assume  that 


f (x  )  =  9  f(x  )  =  0 


F  (!_  £)  =  ik  F(f), 
3x, 


P(l_2f)  =  -k2  F(f) 


Also  define 


3L(f)  =  e"SL  f(t)  dt 


as  the  Laplace  transform.  The  initial  conditions  are  assumed  to  be. 


-  6 


u  (x  ,x  ,0) 

1  1  3 


3  u  (x.,x  ,0)  =  0 
3t  1  1  3 


(1.5) 


u  (x  ,x  ,0) 

3  1  3 


3  u  (x  ,x  ,0)  =  0 
3t  3  1  3 


Letting  flF(u  n=  u  (k,x  ,s),T[f(u  )]=  u  (k,x  ,s),  we  have  the  following 
113^3  33 

transformed  equation: 


C  ( -k 2 )  u  +  C  u  +  (C  +C  )  (ik  )  u  =  +s2  u 

111  1  551,33  1355  1  3,3  1 


C  u  +  C  (-k2)  u  +  (C  +C  )  (ik  )  u  =  +s2  u 

333,33  55  1  3  1355  1  1,3  3 


and  at  x  =  0  ,  the  transformed  boundary  conditions  become 

3 


(1.6) 


u  +  ik  u  =  0 
1,3  13 


(1.7) 


C  u  +  C  ik  u  =  -p  f(k  )g(s) 

333,3  13  11  O  1 

Since  we  are  expecting  surface  wave  propagation,  we  seek  the  solution  of 
(1.6),  (1.7)  in  the  forms: 


u 

1 

il 

M 

i 

u 

3 

- 1 

CM 

-0- 

_ L 

-p(k,s)x 
c  3 


(real  (p)  0) 


(1.8) 
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Therefore  *  the  equations  for  <f>^,  $  are: 


r  .2 


s -  -  C  k2  +  P2  C  -ik  p  (C  +C  ) 

111  55  I  1355 


•ik  p(C  +C  ) 

1  r  1  3  55 


-  s2  -  k2  C  +  p2  C 
155  : 


t — 

1 

4) 

i 

2 

=  0 


(1.9) 


det 


or  for  a  non-trivial  solution, 

=  C  C  p4  +  [C  (-s2-C  k2)  +  C  (-s2-C  k2)  +  k2 (C  +C  ) 2]  p2 

3355^  33  111  55  551  11355 

+  (s2+C  k2) (s2+C  k2)  =  0  (1.10) 

^  111  551 


We  will  choose  the  p's  with  positive  real  parts  to  ensure  the  decay  in 

direction  of  the  surface  wave.  Let  the  solutions  be  p  =  p  ,p  ,  therefore, 

i  2 

we  have 


.0) 


(l) 


p  t  p  :  (j>  =  C  (k  ,s) ,  (j>  - 

y  F1  1  11  2 


i  [ - s 2 -C  k2  +  C  ^  ] 

_ _ 111  551 

k  p  (C  +  C  ) 

11  13  55 


<f>^  =  Ip  C 

1  3  11 


(1.11) 


(2) 


(2) 


p  =p2:  ^  2  =C!(k]|S),  ^  -  CJ;C2  = 


i[-s2-C  k2  +  C  p2] 

_ _ 111  552 

k  p  (C  +C  ) 

12  13  55 


(2) 


Therefore,  the  displacements  have  the  forms: 


U 

1 

=  C  j( kj ,s) 

1 

e  P  iX  3  +  C  (k  ,s) 

2  1 

1 

U 

p 

m  3^ 

3  1 

3  2 

(1.12) 
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From  the  stress-strain  relations 


t  =  C  (u  .  +u  ),t  =Cu  +C  u  , 

13  551,3  3,1  11  111,1  13  3,3 


t  =  C  u  +  C  u 

33  33  3,3  13  1,1 


(1.16) 


we  obtain  the  transforms  of  the  stresses 


t  =  C  { (-p  +ik  ip  )  C  e  Pi  3 

13  55  1  131  1 


+  (-p  +ik  ip  )  C  e  ^2  3} 

2  1  3  2  2 


t  =  (C  ik  -C  p  ip  )  C  e  Fi  3 

11  11  1  .13  13  1  1 


(1.17) 


+  (C  ik  -p  \p  C  )  C  e‘P2X3 

XI  123213  2 


t  =  Op  ip  C  +ik  C  )  C  e  F1  3 

33  13133  113  1 


+  Op  Ip  C  +C  ik  )  C  e  *2  3  , 

23233  13  1  2 


10 


(1.19) 
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II.  EDGE  IMPACT  OF  PLATE  WITH  EDGE  PROTECTION 

To  prevent  failure  of  composite  fan  blades  under  impact  forces,  leading 
edge  protective  strips  have  been  employed.  In  practice,  these  strips  of 
stainless  steel  are  wrapped  around  the  leading  edge.  To  model  this  device, 
we  consider  a  beam  bonded  to  the  edge  of  an  anisotropic  plate  (Figure  l)* 

The  effect  of  the  beam  will  be  to  thwart  the  force  of  impact,  thereby 
decreasing  the  normal  stresses  in  the  composite.  However,  we  shall  show 
that  with  such  a  reduction  in  normal  stress , sizeable  interface  shear  stresses 
can  be  induced. 

With  the  introduction  of  a  beam  of  thickness  b  on  the  edge  of  the  composite 
plate,  the  Rayleigh  wave  behavior  will  depend  on  the  ratio  of  the  wavelength 
to  thickness  ratio  of  each  Fourier  component  in  the  x^  direction.  Thus  one 
should  expect  the  Rayleigh  wave  speed  to  vary  with  b/a,  the  thickness  to 
impact  footprint  ratio.  In  addition  the  Rayleigh  wave  will  become  distorted 
as  it  propagates. 

To  solve  the  edge  strip  problem  the  solution  in  the  composite  plate 
follows  the  same  procedure  as  the  no- strip  case  except  for  the  boundary 
conditions  on  the  edge.  In  place  of  the  zero  stress  conditions  on  the  edge 
we  relate  the  edge  stresses  t^,  tJ3  to  the  motion  of  the  beam  strip.  If 
one  considers  a  small  element  of  the  beam-strip  along  the  x.^  direction,  the 
momentum  balance  equations  in  the  x^  x^  directions  become,  (for  a  plate 
of  unit  thickness) 

pb  ify.  =  Eb  ifu  +  b13 

3t2  3x2 


(2.1) 
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pb  32W  =  -  El  34W  +  I  b  34W  +  b  3t  _  +  t. 


(2.2) 


3x«t‘ 


In  these  equations  U,  W  are  the  x^,  x^  displacements  of  the  beam  element  at  the 
half  thickness,  and  tgg,  t^3  are  the  interface  stresses. 

We  choose  the  compatibility  conditions  between  the  beam  and  plate  dis¬ 


placements 


W  =  Ug,  on  Xg  =  0. 


(2.3) 


U  =  Uj  +  b  3u3  ,  on  Xj  =  0 


(2.4) 


In  the  above  equations  b  is  the  depth  of  the  strip,  E,  I,  Ip  are 
respectively  the  Young's  modulus,  moment  of  inertia  and  rotary  inertia.  Also 
p  f(t)  g(x1)  is  the  edge  loading  now  applied  to  the  outer  protective  strip 
surface. 

The  equations  for  the  plate  remain  as  in  the  free  edge  case  and  a 
solution  is  obtained  by  taking  a  Laplace  transform  on  time  and  a  Fourier 
transform  on  the  space  variable  x j.  With  nondimensionalization  the  solution 
in  the  plate  is  assumed  in  the  form  of  01.12).  The  transform  of  the  plate 
displacements  are 


e'PlX3  +  C, 


e'P2X3 


(2.5) 


where  p^,  p^  are  defined  in  (1.10)  and  ^31 5  $32  are  &iven  (1«11)- 
C2  are  determined  from  the  edge  boundary  conditions.  However,  in  place  of  the 
free  edge  conditions  (1.2)  we  use  the  equations  of  motion  for  the  strip  (2.1), 


(2.2).  Cj,  C2  are  then  solutions  to  the  algebraic  equations 


where 

G  =  -k2Ew  +  p  C  -  pws2  -  ip  (ik2w2E/ 2  +  ikC  +  pw2s2ik/2) 

l  l  55  31  55 

H  =  -k2Ew  +  p  C  -  pws2  -  ip  (ik3w2E/2  +  ikC  +  pw2s2ik/2) 

l  2  55  32  55 

G  ="ikC  +  ikwp  C  J2  +  ip  (  p  C  +  Ew2k**/12  +  wk2C  /2 

13  ri  55  31  l  33  55 

+  k2s2pw3/|2  +  pws2) 

H  =-ikC  +  ikwp  C  J2  +  ip  (  p  C  +  Ew-2k4/12  +  wk2C  /2 

13  '  5o  32  233  55 

+  k2s2pw3/|2  +  pws2) 

where  w  =  b/£,  p,  E,  are  nondimensionalized  quantities  and  p^, 
ip  ,  tp  are  defined  in  (1.11). 
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III.  NUMERICAL  INVERSION 

The  inversions  are  accomplished  by  the  Fast  Fourier  Transform  (FFT) 
techniques  [2],  which  consists  of  a  transformation  from  Laplace  to  Fourier 
transforms,  and  a  two-dimensional  numerical  inversion  using  the  usual 
FFT  alogarithm.  Notice  the  Laplace  inversion  formula 

•  C+i00 

fw  =  ^  f(s>  «st  ^ 

<*  C-i00 

r 

Set  s  =  C  +  ia 


f(t)  =  Lf  f (C+ia)  eCt  eiat  da  (3.1) 

J  _co 

r 

where  C  and  a  are  both  real,  and  C  is  greater  than  the  largest  real 
part  of  all  singularities  of  f(s).  Numerically,  the  double  Fourier  transform 
(or  inversion)  has  the  following  form  [1]: 


f(x,t) 


-i[KxCl-|u  *  Kt(l-±)t] 


N  M 

£  £  f(I,J)  e 

1=1  J=1 


(I-l)x  (J-l) 
27Tl[_N -  +  M 


(3.2) 
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where  N,  M  are  the  number  of  points  in  x  and  t  direction  respectively, 
and  K  ,  K  are  the  corresponding  half-frequency  range. 

X  "C 

For  the  present  problem,  the  determination  of  C  is  through  the  following 
considerations : 

The  form  of  inversion  integrals  are,  in  general. 


p  =  p(k^  ,s)  ,  Re(p)jK) 

(3.3) 

and  A  is  given  by  (1.14).  It  is  easy  to  see  the  singularities  of  the  integral 
of  (2.3)  are: 

a)  Poles  at  A  =  0, 


27Ti 


f  L+  loo 


C-i00 


FCk^s) 
A(k^ ,s) 


-px 


3  St 

e  as 


A  =  (p  -p  )  ik  (C  -ip  ip  C  )  +  (ip  -ip  )  (k2C  +p  p  C  ) 

2  1  1  13  31  32  33  32  31  1  13  ^1  2  33 


(  a  4) 


implies 


(p  -p  ) [C  k2p  p  (C  +C  )2  +  C  (C  p2-s2-C  k2) (C  p2-s2-C  k2)] 

^2  ^1  L  13  1  Ir2  13  55  33  55  1  11  1  55  2  11  1 


(C  +C  )[p  (C  p2-s2-C  k2)  -  p  (C  p2-s2-C  k2) ] (k2C  +  p  p  C  ) 

1  3  5  5  1  5  5*  2  1  1  1  r2  55ri  11  1  1  13  r  lr 2  33 


0 


(3.5) 
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Interpretation  of  this  condition  is  best  understood 

for  the  particular  case  of  an  isotropic  material,  i.e.  for 


C  =  A  +  2y  =  C  ,  C  =  A,  C  =  y, 

11  33  13  55 

choose  C  =  C  =  y 

6  6  5  5 

is  easy  to  see  that  from  (1.10),  (1.11) 


det  =  [yp2  -  (yk2  +  ps2)]  {(2y+A)  p2  -  [(A+2y)  k2  +  ps2]} 

p2  =  (yk2  +  ps2)/y,  p2  =  [(A+2y)  k2  +  ps2]/(A+2y) 

11  2 

ip  =  ik  /p  ,  ip  =  -p  /ik  , 

31  1*1  32  r  2  1 

Here  [y/p]^2  =  v  ,  is  the  shear  wave  speed  and  [  (A+2y)/p] 1//2  v  ,  is  the 
s  P 

longitudinal  or  pressure  wave  speed  for  isotropic  materials. 


* 


17 


Thus  for  the  isotropic  case  that,  A  -  0  implies 


4yp  +  (p 
2  1 


(A  +  2p)  ]  =  0 


C3.6) 


where  p^,  p^  are  defined  above.  This  is  the  equation  for  the  Rayleigh  wave 

1/2 

speed  v^.=  (is  /k)  which  is  found  as  a  real  root  of  (3.6)  (see  e.g.[3]). 

For  the  case  A  =  y (Poissons  ratio  =  0.25),  vD  =  0.919  v  . 

R  s 

For  the  anisotropic  case  a  computational  scheme  to  calculate  the  zeroes 
of  (3.5)  has  been  written,  (Figure  2). 


b)  Branch  points:  The  branch  points  of  the  integrand  ih(3.3)  Ate  the  same 


as  those  of  the  functions  p  (k  ,s),  p  (k  ,s),  they  are: 

11  2  1 

i)  p  =  0,  or  p  =0  which  implies  (C  k2  +  ps2)  (C  k2  +  ps2)  =  0 
1  2  11  1  55 

i.e.,  k/s  =  ±i/p/C  ,  ±i/p/C  pure  imaginary. 


55 


1 1 


These  correspond  to  longitudinal  and  shear  wave  speeds  for  an  isotropic  material 


ii)  p 


=  + 


/ 


1  ,2 


B  ±  /  B2  -  4AC  / 2A 


with  A  =  c  C  >  0 

3  3  5  5 


(3.7) 
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B  = 

C  = 


-ps2 (C  +C  ) 

3  3  5  5 


+  k2 [ (C  +C 

1  3 


e  e  -c2  ] 

11  33  55 


(C  k2+  ps2)  (C  k2+Ps2) 

11  5  5 


These  branch  points  are  those  values  of  s/k  which  render  B  -  4AC  =  0, 
and  are  branch  points  of  second  order.  ^  The  distribution  of  these  points  is 
shown  in  Figure  3.  It  has  been  shown  [5]  that  the  contribution  of  these  branch 
points  to  the  value  of  the  integrals  (3.3)  is  important  only  when  one  considers 
the  multi-reflected  and  refracted  waves  in  layered  media,  or  when  the  position 
of  interest  is  very  close  to  the  impact  origin.  In  this  study  we  were  more 
concerned  with  how  the  energy  is  propagated  away  from  the  impact  point,  which 
is  mainly  associated  with  surface  waves,  thus  we  ignored  the  contribution  of 
these  branch  points. ^ 

The  contour  of  integration  for  the  Laplace  inversion  is  as  shown  in 
Figure  3.  Notice  the  branch  cuts  are  extended  to  negative  infinity,  in 
accordance  with  the  requirement  that  C  >  max  real  part  of  the  singularities. 
The  requirement  that  real  (p)  0  also  determines  the  correct  sheet 

of  the  Riemann  surfaces.  Numerically,  since  the  branch  points  are  located 

ct 

at  s/k  =  constant,  as  k  gets  large  C  should  be  large,  and  the  factor  e 
in  the  Laplace  inversion  expression  will  rise  sharply  to  an  unmanageable  size. 
Since,  in  the  last  paragraph,  we  have  noticed  the  contribution  of  the  branch 
points  is  unimportant, ' a  path  r2  is  chosen  to  replace  Tj  by  the  Cauchy's 
integral  theorem.  Notice  the  advantage  of  integration  along  r2  is  that  Cq 
is  a  positive  constant  independent  of  k.  The  determination  of  optimal  Co 
is  discussed  in  [2].  Here,  in  order  to  minimize  the  aliasing  and  round  off 
errors  in  numerical  computations  simultaneously,  we  choose 
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IV.  RESULTS 

A  computer  program  has  been  written  to  calculate  the  stresses  in  a 
plate  with  an  elastic  beam  on  one  edge  under  a  transient  impact  load 
distribution  along  the  edge.  A  program  description,  flow  charts,  input 
data  formats,  and  sample  printout  of  the  program  are  contained  in  the 
appendices  to  this  report.  In  this  section  we  will  summarize  some  of 
the  results  obtained  from  this  computer  program.  These  results  were 
calculated  for  an  anisotropic  plate  with  effective  elastic  constants 
of  55%  graphite  fiber/epoxy  matrix  composite  obtained  from  Reference  7 
and  summarized  in  the  Table. 

No- strip  case: 

The  Rayleigh  wave  can  be  seen  in  the  stress  t^  on  the  edge  as  shown 
in  Figures  4,  5  for  a  graphite  fiber-epoxy  composite  for  layup  angles  0,  _+ 
15°.  After  the  initial  contact  time,  the  stress  is  observed  to  propagate 
with  little  change  at  a  speed  near  the  calculated  Rayleigh  speed 
(Figure  2).  This  wave  can  be  observed  in  the  computed  output  in  the 
space- time  (x^,  -t)  plane  Figure  7,  as  a  band  of  non-zero  values  along 
a  diagonal  from  the  upper  left  to  the  lower  right  comer  of  the  Xp  -t 
plane.  Caution  is  urged  in  using  this  program  since  spurious  waves 
can  enter  the  calculations  due  to  the  periodic  nature  of  the  finite 
numerical  Fourier  Transform.  These  spurious  waves  are  data  bands  which 
lie  along  diagonals  from  upper  right  to  lower  left.  In  other  words, 
only  disturbances  emanating  from  the  impact  source  in  the  upper  left 
corner  of  the  x^  -t  plane  of  Figure  7  should  be  valid.  A  computer  map 
of  the  space  time  history  of  the  edge  impact  stress  is  shown  in  Figure  6. 
A  contact  time  of  35  ysec  and  contact  length  2  cm  was  used  in  these  cal¬ 


culations  . 
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As  is  characteristic  of  surface  wave  effects,  the  stresses  due  to 

impulsive  loading  on  the  edge  decrease  with  distance  from  the  edge.  This 

is  shown  in  Figure  8  for  two  different  layup  angles.  The  normal  stress 

t_,  appears  to  decrease  to  about  1/4  of  its  value  on  the  surface  at  a 
o  o 

depth  equal  to  one  half  of  the  loading  length  a.  The  rate  of  decay  from 
the  edge  depends  on  the  layup  angle.  Another  characteristic  of  edge 
impact  is  the  development  of  tension  in  the  normal  stress  t33  under  the  impact 
point.  This  is  shown  in  Figure  9.  Thus  while  the  compression  part  decreases 
with  distance  from  the  edge  a  tension  tail  developes  in  the  wave. 


The  effect  of  layup  angle  on  the  impact  stresses  can  be  seen  in 
Figures  10,  11.  The  stress  t]x  at  the  edge  is  larger  than  the  impact 
pressure  and  decreases  as  the  layup  angle  goes  from  0°  to  _+  45°  (Figure  10). 

Below  the  surface  or  edge,  the  peak  normal  stress  t^  at  x^  =  0  is 

a  minimum  for  layup  angles  near  _+  30°,  while  the  shear  stress  t^3  increases 
as  the  fiber  angle  goes  from  0°  to  +_  45°  (Figure  10) . 

An  unexpected  result  is  the  shift  of  the  maximum  normal  stress  t33  to 

points  off  the  impact  axis  =  0  for  layup  angles  greater  than  about 
_+  30.  A  pronounced  peak  in  t33  versus  x^  beyond  the  impact  pressure 

foot  print,  can  be  seen' in  Figure  11  for  _+  45°  layup  angle. 

Impact  protection  strip  case: 

The  effect  of  bonding  an  edge  impact  protection  strip  to  the  half 
plane  is  shown  in  Figures  12-16.  In  Figures  12,  13,  the  increase  in  the 
thickness  of  a  steel  strip  produces  a  decrease  in  the  interface  stresses 
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tjj,  but  creates  an  interface  shear  stress  at  the  strip-composite 

interface.  This  shear  stress  reaches  a  maximum  for  strip  thickness  less 
than  the  impact  footprint  length  and  decreases  for  greater  strip  thick¬ 
nesses.  Thus,  if  the  strip  is  too  thin,  debonding  can  occur  under  im¬ 
pact  due  to  induced  interface  shear. 

In  Figure  14  one  can  see  that  increasing  the  strip  thickness  decreases 
the  peak  normal  stress  t^g  and  redistributes  the  load  over  a  longer 
length  under  the  strip.  However,  while  the  peak  compression  stress  is 
decreased  by  the  strip,  tension  is  created  which  could  also  produce 
debonding  of  the  strip  from  the  composite. 

For  the  no- strip  case  a  Rayleigh  wave  was  seen  to  propagate  along 
the  edge  relatively  unperturbed  (Figures  4,  5).  With  the  strip  present, 
(Fig. 15)  this  wave  becomes  distorted  as  time  increases.  In  fact,  the 
beam- strip  boundary  conditions  introduce  dispersion  in  the  edge  waves 
which  make  the  Rayleigh  edge  wave  velocity  dependent  in  the  effective 
wave  length  of  the  disturbance. 

Finally  in  Figure  16  shear  stress  distributions  along  the  strip- 
composite  interface  are  shown  for  a  thickness  near  the  shear  peak 
(b/a  =  0.25)  and  another  for  b/a  =2.0.  In  the  latter  case  the  shear 
is  distributed  over  a  larger  length  resulting  in  a  lower  peak  stress. 

Also  the  strip  delays  the  time  of  maximum  shear  from  1/2  to  3/4  Tg 
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Summary  of  Results 

An  analytical-numerical  method  has  been  developed  to  solve  the  response 

of  a  composite  plate  with  a  bonded  edge  strip  to  in-plane  impact  type 

forces  on  the  edge.  Results  of  computer  simulations  reveal  the  following: 

1)  Rayleigh  edge  waves  can  propagate  away  from  the  impact  site  with  tension 

and  compression  up  to  values  of  the  impact  pressure,  depending  on  layup  angle. 

2)  Normal  to  the  edge,  the  initial  peak  compression  pulse  decreases  as 
it  propagates  into  the  plate  but  a  tension  tail  develops  as  it 
propagates  away  from  the  impact  site. 

3)  The  edge  stress  t^  under  impact  is  decreased  as  the  fiber  layup  angle 

goes  from  0°  to  +_  45°. 

4)  Protection  strips  of  thickness  less  than  half  the  impact  length  can 
develop  large  interface  shear  under  impact. 

5)  The  normal  and  edge  stresses  t33,  tn  at  the  edge  can  be  decreased 
significantly  by  protection  strips  of  thickness  greater  than  the  half 
impact  length. 
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APPENDIX  A 

The  Determination  of  Parameter  C 

0 

Consider  the  Laplace  inversion  of  a  function  g(t)  as 

(•Co+i°° 

g(t)  =  2ii  gO)  eSt  ds 

Co‘io° 

and  let  s  =  C  +iu),  C  and  to  are  real.  We  can  change  (A.l)  into  the 
o  o 

Fourier  inversion  formula 

C  t  f  °°  .  C  t 

g(t)  =  ~  e  °  g (Co+itl»)  e1Wt  dua  =  e  0  x(t)  (a.  2) 

-OO 

which  is  then  inverted  by  the  Fast  Fourier  Transform  technique.  In  the 
FFT  scheme,  a  continuous  function  g(CQ+ito)  is  discretized  and  the 
infinite  interval  of  integration  is  truncated.  The  error  due  to  trunca¬ 
tion  depends  on  each  problem  but  doesn't  depend  on  Cq,  thus  in  the 

determination  of  C  ,  we  will  assume  the  truncation  error  is  negligible. 

o 

It  is  shown  [2]  that  the  discretizing  of  the  transform  in  one  domain  will 
cause  aliasing  error  in  the  other  domain,  e.g.  sampling  g  at  N  points  in  a 

frequency  interval,  0  <  »  <  12  ,  will  produce  a  transformed  function  x^ft) 
which  is  periodic  and  which  differs  substantially  from  x(t)  for  large 
enough  t.  For  even  x(t)  this  difference  can  be  shown  [2]  to  be  given  by, 

x  (t)  =  x(t)  +  x(t-T) 

for  0  <  t  <  T/2  where  T  =  N/ft. 
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It  has  been  shown  [2]  also  that  the  aliasing  error  is  approximated 
,C  (T-2t) 

by  E  (t)  =  e  °  g(T-t)  for  the  Laplace  inversion  (A. 2). 


Notice 


the  aliasing  error  is  a  decreasing  function  of  C 


The  other  source  of  error  is  of  course  the  round  off  error  in  compu- 

V 

tation.  Since  we  multiply  the  resulting  x(t)  by  e  to  get  g(t), 
the  rounding  error  is  of  the  form 


C  T 

Er (t)  =  e  r (t) . 


The  error  bounds  are  then 


e  =  I  Max  Ea(t) |  =  e 


-C  (T-2t) 

MaxJg(T-t) 


e  =  | Max  Er(t)|  =  e  Max|r(t) 

2 


Equating  and  e^,  the  optimal  Cq  is  then 


Jin  (Max  g(T-t)/Max  r(t)} 


Chosing  T  =  T/4,  therefore 


0  <  t  <  T 


C  =  w  &n(g/r)  , 
o  5 1 


Max  g(T-t) ,  r  =  Max  r(t) 


Notice,  empirically,  f  =2  ^  10"6  on  single  precision  IBM  360 


systems , 


APPENDIX  B 


A  FLOW  CHART  OF  THE  PROGRAM 
1.  Main  Program 


Read  in  material  data, 
strip  data 


Calculate 

speeds 

wave 

\ 

f 

|  Calculate 

1  T 

max 

X  and 

max 

|  Non-dimensionalization 


Calculate  Laplace 
Inversion  Parameter 


Calculate  singularities 
including  Rayleigh  speed 


Calculate  the  trans¬ 
formed  expressions 


Plot  the  Relative 
Values 


Repeat  for  different  x. 
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APPENDIX  C 

NOTES  ON  COMPUTER  PROGRAM 

The  choice  of  scales  is  very  essential  to  the  success  of  the  present 
computational  method.  It  is  noticed  that  the  accuracy  depends  on  the 
number  of  points  employed  and  the  range  of  frequency  spectra  covered. 
Considering  limitations  of  both  computer  storage  and  time,  a  time-space 
grid  of  32  x  64  points  was  chosen  for  t  >  0,  x  >  0.  Thus,  the  non- 
dimensionalization  of  all  equations  and  quantities  are  both  necessary  and 
important  to  the  obtaining  of  meaningful  data  from  the  limited  grid  size. 
The  numerical  inaccuracies  introduced  have  several  origins: 

1.  Theoretically,  error  has  been  introduced  by  the  neglecting  of  the 

outer  branch  points  contributions.  It  has  been  shown,  for  isotropic 

-  2 

cases,  the  contributions  of  these  branch  points  behaved  like  r  at 

large  distance  from  a  delta  function  loading  at  the  origin  r  =  0. 

1 

Compared  with  the  r  2  decreases  of  the  contribution  of  the  residue. 

Thus,  for  small  r,  or  near  the  origin,  the  errors  might  be  significant. 

An  asymptotic  form  of  the  behavior  of  small  r  has  been  deduced  for  a 
simple  delta  loading  at  origin  on  an  isotropic  half-space  [  ] .  It  is 

shown  the  error  thus  introduced  is  of  the  order  of  5%  maximum  response 
in  stress. 

2.  The  aliasing  error  introduced  through  the  periodizing  of  the  functions. 

3.  The  round-off  error  in  Laplace  inversion  along  with  aliasing  error . 
have  been  discussed  in  the  determination  of  C^  (Appendix  A) •  It  is 
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found  that  error  3  is  more  serious  of  the  two.  Hence  in  calculations  the 
maximum  non-dimensionalized  time  should  be  restricted  to  below  6  or  8 
for  reasonably  good  results. 

4.  Errors  due  to  reflections  at  the  boundary.  Since  the  space  grid  is 
finite,  it  has  been  observed  that  whenever  a  wave  hit  the  boundary  of 
chosen  space,  a  sizable  numerical  error  will  start  propagating  in,  as  if 
the  wave  were  reflected  from  the  boundary.  The  basic  reason  is  due  to 
the  periodization  of  the  space  (x^)  domain.  In  computation,  this 
error  should  be  avoided.  To  correctly  determine  the  extent  of  the  space 

(x  )  domain,  a  priori  recognization  of  significant  wave  speed  (at  which 

X 

most  of  the  energy  travels)  is  important.  Usually  an  estimation  will  be 
sufficient.  Then  the  proper  nondimensionalizing  constants  can  be 


chosen. 
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APPENDIX  D 

PROGRAM  INPUT  AND  OUTPUT 

A.  Method:, 

FFT  alogarithm  for  numerical  inversion 


B.  Input  Data  Cards: 

Card  1.  NTEST,  NSTRIP,  NP  (315) 


NTEST  = 

1 

A  test  program 

for  FFT  (2-D) 

0 

The  present 

program 

NSTRIP  = 

1 

With  strip 

0 

No  strip 

NP 

1 

Calculate : 

u 

i 

only 

2 

V 

u 

3 

3 

U  , 
1 

U  ,  t 

3  3  3 

4 

V 

u  ,  t  ,,t 

3  3  3  1  1 

5 

V 

u  ,  t  ,  t 

3  3  3  1  1 

Card  2.  CC(I)  I  =  1,9,  RHO,  ANGLE  (8E10,4) 

CC(I=1,9)  Material  constants  of  composite,  in  the  order 

C  ,  C  ,  C  ,C  ,C  ,  C  ,  C  ,  C  ,C  .  (psi) 

11  22  33  kk  55  66  12  13  23 

RHO  Density  of  composite  (fe/cm3) 

ANGLE  (degrees)  lay-up  angle 

Card  3.  VEL  ,  DM,  El,  ANU,  DEN  (8E10.4) 

VEL  -  Velocity  of  incoming  particle  m/sec 

DM  -  Diameter  of  impacting  object  cm 

El  -  Youngs  Modulus  of  impacting  object  (psi) 

ANU  -  Poissons  ratio  of  impacting  object 

DEN  -  Density  of  impacting  object  (gm/cc) 
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Card  4.  NSTRES ,  NK3,  DX3  (215,  F10,3) 

NSTRES  =  1 

NK3:  Number  of  steps  in  direction 

DX3‘  Ax  /£,  step  size  in  x  direction 

3  3 

Photo  elastic  fringe  order  computer  map  in  the  x-t  plane 


Card  5. 


RO,  W,  ES 

(3E10,4) 

RO: 

density  of  strip  beam  (g/cm3) 

W: 

depth  of  strip  beam  (cm) 

ES:. 

Young's  modulus  of  strip  beam  (psi) 

C.  Output: 

1.  Test  problem:  Appendis  III 

2.  Values  of  displacements  1  :  u^ 

2  :  u 

3 

(stress)  3  :  t 

4  :  t 

1  3 

5  :  t 

i  i 


3.  Relative  magnitude  computer  maps  of  displacements  and  stresses 


- ^  x  /£ 

1 

si' 

t/t* 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


********  program  to  calculate  stresses  due  to  edge  impact  of  a  plate* 

THIS  PROGRAM  CALCULATES  THE  ELASTIC  RESPONSE  OF  AN 
ANISOTROPIC  PLATE  TO  AN  IN-PLANE  EDGE  IMPACT  FORCE  ON  X3=0.0 
WHEN  NSTRIP  =1  ,  THE  PROGRAM  PLACES  AN  IMPACT  PROTECTION  STRIP 
,CR  ELASTIC  BEAM  ON  THE  EDGE.  THE  IMPACT  FORCE  IS  A  HALF  SINE 
FUNCTION  IN  TIME  AND  IS  NON  ZERO  FOR  0<T<T0# (MICROSEC) .  THE 
FORCE  IS  DISTRIBUTED  ALONG  THE  EDGE  AS  PQ(1-X1**2),  WHERE 
XI  IS  NORMALIZED  BY  THE  HALF  WIDTH  OF  THE  IMPACT  CONTACT 
LENGTH. 

THE  METHOD  EMPLOYS  A  FOURIER  TRANSFORM  IN  THE  EDGE  DIRECTION 
X 1 ,  AND  A  LAPLACE  TRANSFORM  IN  THE  TIME  DIMENSION.  THE  TRANSFORM 
OF  THE  FORCING  FUNCTION  IS  GIVEN  IN  THE  SUBROUTINE  CALC UL  AND 
THE  SOLUTION  IS  OBTAINED  USING  A  2-DIMENSIONAL  FAST  FOURIER 
INVERSION  ROUTINE  CALLED  'FOURT' • 

THE  OUTPUT  FOR  A  GIVEN  DEPTH  X3  CONSISTS  OF  DISPLACEMENTS 
D  1 , U3, AND  STRESSES  T33,T11,T13,  IN  THE  Xl-TIME  PLANE  . 

THE  STRESSES  ARE  NORMALIZED  BY  THE  ELASTIC  CONSTANT  C(6,6). 

IMPACT  FORCE 
I 
I 
I 
I 

V  I  V 
V  I  V 
VI  V 
V 


STRIP  »>  ***************  «<  STRIP 


C  C - >  XI 


C 

c  I 

C 

COMPOSITE 

C  l 

c 

C  I 

c 

HALF  SPACE 

C  l 

c 

C  V 

c 

C  X3 

INPUT  DATA 

NN  ( 1) — 2.0*  MAX  XI  DISTANCE,  NN(2)--MAX  NO*  OF  TIME  UNITS 
NSTRIP=0,NO  STRIP, .. NSTR I P=1 , WITH  STRIP 
NP=1,--,5,  CALCULATES  U1 , 03, T33 ,T1 1 ,T1 3 , IN  THAT  ORDER 
NP=6,  CALCULATES  DISPL.  OF  A  BEAM  ON  AN  ELASTIC  FOUNDATION 
CC  (9)  ,  ELASTIC  CONSTANTS  OF  PLATE  IN  THE  ORDER 

C11 ,C22,C33,C44,C55,C66,C23,C13,C12, IN  PS! 
RHO, DENSITY  OF  PLATE  IN  UNITS  GM/CC 

ANGLE--LAYUP  ANGLE  OF  COMPOSITE  PLATE, DEG . /FOR  INFO  ONLY 
VEL,  VELOCITY  OF  INCOMING  OBJECT  METERS/SEC 
DM,  DIAMETER  OF  IMPACTING  OBJECT  -CM. 

E 1 , ANU ,  YOUNG'S  MODULUS  AND  POISSON'S  RATIO  FOR  IMPACTING  BODY 

DEN,  DENSITY  OF  IMPACTING  BODY 

NSTRESS=1 

NK3,  NO . OF  DEPTHS  X3,  (FOR  NK3^1,X3=0) 

DX3,  INCREMENT  IN  DEPTH  X 3 (NORMALIZED  BY  AO) 

RO,  DENSITY  OF  PROTECTIVE  STRIP  GM/CC 


*  This  program  has  two  extra  cards  to  override  the  Hertz  contact  time  and 
contact  length  calculation.  Remove  cards  #51,  52. 
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ooei 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 


C 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


w ,  THICKNESS  TO  WIDTH  RATIO  OF  BEAM 
ES,  YOUNG'S  MODULUS  FOE  BEAM 

WL — FOURIER  WAVELENGTH  (CM)  OF  THE  ORDER  OF  AO  OR  LESS 
WT — FOURIER  WAVE  PERIOD  (SEC)  OF  THE  ORDER  OF  TO  OR  LESS 
CHOICE  OF  ML, WT  DETERMINES  DX  f  DT — DX=WL/2, DT~ WT/2 
FG#  TRANSFORM  OF  NORMALIZED  FORCING  FUNCTION  F(X1)*G{T) 

THIS  IS  PROVIDED  IN  PROGRAM  BUT  CAN  BE  CHANGED  BY  THE  USER 

OUTPUT  DATA 

TC,TC  CONTACT  TIME  (SEC,1.E~6  SEC  )  (FROM  HERTZ  THEORY) 

AO, A,  HALF  THE  IMPACT  CONTACT  LENGTH  ,CM 

FO,  MAX  IMPACT  FORCE  FROM  HERTZ  THEORY  , NEWTONS 

E ,  WAVE  SPEED  IN  PLATE  SQRT  (C66/RHO) 

OR  WAVE  SPEED  IND  BEAM  STHI P, SQET ( £1 /DEN)  UNITS  CM/SEC 
CL-SQRT  (Cl 1/RHO)  LONGITUDINAL  WAVE  SPEED  ALONG  EDGE  IN  PLATE 
CS^SQBT  (C55/RHO)  ,  SHEAR  SPEED  ALONG  EDGE  OF  PLATE f CM/SEC 
CR,  RAYLEIGH  WAVE  SPEED  ALONG  FREE  EDGE  OF  PLATE 

DX#  DT  SPACE  TIME  INCREMENTS  IN  Xl-T  SPACE  UNITS — CM  AND  SEC 
DATA  (I, J)  FORMALIZED  TRANSFORM  OF  ONE  OF  DISPL.  OR  STRESSES 
•-BEFORE  CALL  FOURT,  AFTER  CALL  FOURT  DATA  IS  A  2  DIM  MATRIX  OF 
DISPL.  OR  STRESSES  IN  Xl-T  SPACE  f DEPENDING  ON  VALUE  OF  K 
IN  THE  LOOP  1  DO  4  K=1,NPf 

Ul,U3, NORMALIZED  DISPL.  IN  PLANE  OF  PLATE  (E.G.  Ul/AO) 
T33#*T11,T13,  NORM.  STRESSES  (E.G.  T33/C66) 

— NOTE--  T33  ON  X3=0  SHOULD  REPRODUCE  THE  FORCING  FUNCTION 
WHEN  THERE  IS  NO  STRIP 

--NOTE--  AS  A  CHECK  T13=0  ON  X3=0  WHEN  THERE  IS  NO  STRIP 
THE  FPINGE  ORDER  MAP  PLOTS  THE  DIFF  IN  PRINCIPAL  STRESSES  AND 
MAY  BE  USED  TO  COMPARE  WITH  PHOTOELASTIC  EXPERIMENTS  OR  TO 
LOOK  FOR  POINTS  OF  MAX  IN  PLANE  SHEAR  STRESSES 

**********  ********** 

THIS  PROGRAM  HAS  BEEN  WRITTEN  BY  F.MOON  AND  C-K  KANG  UNDER 

A  GRANT  TO  PRINCETON  UNIVERSITY  FROM  THE  NASA  LEWIS  RESEARCH 
LAB. 

********** ********** 

--NOTE  TO  THE  USER —  IN  THE  OUTPUT  MAPS  OF  STRESSES  OR  DISPL., 
YOU  WILL  NOTICE  EANDS  OF  SIMILAR  NUMBERS  RUNNING  FROM  THE 
UPPER  LEFT  CORNER  TO  THE  LOWER  RIGHTCORNER  —  THESE  ABE  WAVES 
WHICH  EMINATE  FROM  THE  IMPACT  POINT — HOWEVER-  WAVES  RUNNING 
FROM  RIGHT  UPPER  TO  LEFT  LOWER  CORNER  ARE  SPURIOUS  DUE  TO  THE 
DESCRETENESS  OF  THE  NUMERICAL  FOURIER  INVERSION  PROGRAM 
— ALSO  DATA  FOR  TIMES  NEAR  TM AX  AT  THE  BOTTOM  OF  THE  MAPS  ARE 
USUALLY  SPURIOUS  AND  SHOULD  NOT  BE  USED 

COMMON  /MMC/  D, DKl 

COMMON  /MC/DK2, FG1 ,T0,Cl 1 ,C 1 3 ,C 33 ,C 55 , RHO , RO, W , ES , NS TRI P, AO 
DIMENSION  DATA  (128,32)  , MM  (40)  ,CC  (9) 

DIMENSION  NN (2) 

DIMENSION  RDATA  (40) 

DIMENSION  C I D  AT  A  (40) 

DIMENSION  FRNGE (64,32) 

DIMENSION  T1 1 (64, 32) ,T33  (64,32)  ,T13  (64, 32) 

COMPLEX  DATA , S 

COMPLEX  Pi, P2,S1,S2,C1,C2,B,C,D,SI 
COMPLEX  DK2,FG, SLAP, DKl 
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0012 

CCMPLEX+16  BS 

0013 

NN  ( 1)  =  128 

0014 

NN (2) =32 

C*** 

READ  IN  AND  WRITE  OUT  DATA  AND  PARAMETERS 

c*** 

RHO,AO,TO  MUST  BE  IN  C.G.S.  UNITS,  PO  MUST  BE  CONSISTENT  HI1 

0015 

CALL  INDUMP 

0016 

READ  (5,102)  NTE ST ,N STRIP , NP 

0017 

.READ  (5,100)  CC,RHO, ANGLE 

0018 

WRITE  (6,490)  CC 

0019 

WRITE (6,491)  BHO, ANGLE 

0020 

READ (5, 100)  V EL, DM, El ,ANU,DEN 

0021 

P0=CC (6) 

0022 

C 1 1=  CC (1) -CC (7) **2/CC  (2) 

0023 

C33=  CC (3) -CC (9) **2/CC (2) 

0024 

C  13=  CC  (8)  -CC  (7)  *CC  (9)  /CC  (2) 

0025 

C55=  CC  (5) 

0026 

PI=  3.14159265 

0027 

RHO=  RHO/6. 895* 1. E- 4 

0028 

READ  (5,101)  NSTRES,NK3,DX3 

0029 

WRITE  (6,504)  NSTRES,NK3,DX3 

0030 

SI=  CMPLX (0., 1.) 

c********  CALCULATE  THE  IMPACT  CONTACT  TIME  , RADIUS, AND  PRESSURE 

C********  BASED  ON  HERTZ  CONTACT  THEORY 

0031 

R=DM/2. 0 

0032 

R=R* 1 . E-2 

0033 

E1=E1*6895. 0 

0034 

E2=CC  (2)  *6895.0 

0035 

DEN=DEN*1 . 0E3 

0036 

AMASS=4./3. *PI*R**3*DEN 

0037 

AK2=4./3.*SQRT(R)*El/( ( 1 . 0-A NU* ANU) +E1/E2) 

0038 

ALF=5./4.*AMASS  *VEL*VEL/AK2 

0039 

ALF=  (ALE) **0.4 

0040 

TC=2 . 943*A  LF/VEL 

0041 

F0= 1 . 14*AM ASS*VEL*VEL/ALF 

0042 

A=SQRT (ALF*R) 

0043 

A= 1 . 0E2* A 

0044 

TC=TC* 1 . 0  E6 

0045 

WRITE  (6, 710)  VEL,TC, A,F0 

C************** ************************ 

c*** 

CALCULATE  THE  NON-DIMENSIONAL  PARAMETERS 

0046 

EE=  CC  (6) 

0047 

E=  SQRT (EE/RHO) 

0048 

RE=  RHO 

c 

DEFINE  TRANSFORM  SPACE  AND  DISTANCE-TIME  SPACE 

c*** 

UNIT  DISTANCE  — CM. 

0049 

A  0=A 

00  50 

T0=TC*1 . E-6 

c 

TEST  CASE  T0==  35E-6,  A0=1  CM 

0051 

A0=  1 . 0 

0052 

T0=35. E-6 

0053 

AL=AC 

c*** 

UNIT  TIME — SEC. 

0054 

TE=AI/E 

c*** 

SMALLEST  WAVELENGTHS 

0055 

W  L=  A  0/ 1 . 5 

0056 

WT=TC/10. 0 
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0057 

0053 

0059 

0060 

0061 

0062 

0063 


0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0C78 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 


0099 
0100 
010  1 
0102 
0103 
0104 
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C***  DIMENSIONS  OF  TRANSFORM  SPACE 

XK=2 . 0*PI*AL/WL  ; 

TK=2. 0*PI*TE/WT 

C***  NCNDIMENSIONALIZE  CONSTANTS 
A0=A0/AL 
T  0=T  0/T.E 
DT=HT/2.0 
DX=WL/2.0 

WRITE  (6,621)  WT,WL,TK,XK 

C********* ***************** ************ 

C  CALCULATE  RAYLEIGH  WAVE  SPEED 

X1=0.0 
X2=C55 
N=0 

DO  90  1=1,100 
N  =  N+ 1 

X= (Xl+X2)/2.0 
D 1=  (C55-X) /  {Cl  1  -X) 

D1=SQRT(D1) 

D2=  (Cl  3)  **2-C33*  (C11-X) 

D3=C55*C33 
D3=SQRT  (D3) 

F=D 1 *D2/D3 
F=X+F 

IF (F)  81,82,83 

81  X1=X 

GO  TO  92 

82  GO  TO  91 

83  X2=  X 

92  CONTINOE 
F1=F/C55 

DIFF=ABS (FI) -1 . OE-4 
IF(DIFF)  91,91,90 

90  CONTINUE 

91  CONTINUE 

WRITE  (6,701)  C11,C33,C55,C1.3 

701  FORMAT  (//,10X,6H  Cll  =,E12.4,6H  C33  =,E12.4,6H  C55  =,E12.4, 

16H  C 13  =  ,E  1 2.  4 ,/) 

CL=Cl 1/RHO 
CL=SQRT (CL) 

CS=C55/RHO 
CS=SQRT  (CS) 

CR=X/RHO 

CR=SCRT(CP) 

RDS=CR/CS 

WRITE (6,702)  CL,CS,CR,RDS,N 

702  FORMAT (//, 10X, 1 2HLONG. SPEED  =,E12.4, 13H SHEAR  SPEED  =,E12.4 
1  // , 1  OX, 16H RAYLEIGH  SPEED  =,E12.4,8H  CR/CS  =, F10. 5, 5X,I5) 

C***************** ********** *********** 

C***  CONSTANTS  FOR  STRIP  CASE 
READ  (5,100)  RO , W , ES 
WRITE  (6,521)  RO, W , ES 
RO=  RO/6. 895*1 . E-4 
IF  (RO)  51,51,50 
50  CONTINUE 

E=  SQRT  (ES/RO) 
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0105 

WRITE  (6,510)  E 

0106 

51 

CONTINUE 

0107 

RO=  FO/EE 

0108 

ES=  ES/EE 

0109 

W=  W/AO/AL 

C* *******  ********  **************** ****** 

0110 

205 

CONTINUE 

01  1 1 

C11=  Cl  1/EE 

0112 

Cl  3=  Cl  3/BE 

0113 

C33=  C33/EE 

0114 

C55=  C 5 5/EE 

0115 

P0=  PO/EE 

0116 

RHO=  BHO/RE 

0117 

FG 1=  P0*T0 

0118 

206 

CONTINUE 

0119 

N=  NN  ( 1) 

0120 

M=  NN  (2) 

c*** 

CALCULATE  THE  LAPLACE  INVERSION  PARAMETER 

0121 

WT=«T/TE 

0122 

RH=  P0*WT*1 . E6 

0123 

CH=  2./3./M/WT*  ALOG (RH) 

0124 

CLAP=  CH 

c*** 

CALCULATE  THE  SECOND  ORDER  BRANCH  POINTS 

0125 

A=  RHO*  (C33-C55) **2*R HO 

0  126 

E=  -2.* RHO* ( (C3  3+C55) * (Cl 3*C 1 3+2 . *C1 3*C55 
1*  (Cl 1  +  C55) ) 

0127 

F=  (Cl 3*Cl 3+2 . *C13*C55-C1 1*C33)  **2-4.*Cl1 

0128 

BS=  1.D0*E*E-1 .D0*4.*A*F 

0129  • 

D=  CDSQRT(BS) 

0130 

P1=  .  5/A* (-E  +  D) 

0131 

P2=  F/P 1/A 

0132 

P1=  CSQRT  (PI) 

0133 

P2=  CSQRT (P2) 

0134 

WRITE  (6,509) 

0135 

WRITE  (6,506)  P1,P2 

0136 

204 

WRITE  (6,508)  CLAP 

0137 

K2=  N/2+1 

0138 

M  2=  M/2+1 

0139 

IF  (NTEST.EQ.1)  GO  TO  211 

0140 

211 

IF  (NSTRIP.EQ. 1)  WRITE  (6,514) 

c*** 

GENERATE  THE  TRANSFORMED  EXPRESSIONS 

0141 

REWIND  2 

0142 

DO  1  I-  1 ,  N 

0143 

DR1=  2.*XK/N* (I-.5) -XK 

0144 

DO  1  J=  1,M 

0145 

DK2F=2 . *TK/M* ( J- . 5) -TK 

0146 

SL AP=  CHPLX (CLAP,DK2F) 

0147 

DK2=  -SI*SLAP 

0148 

IF  (NTEST.NE.1)  GO  TO  201 

0149 

FG=  PI*  (1.+CEXP  (-SLAP)  )/(SLAP*SLAP+PI*PI) 
1C SIN (DK1)/DK1* (1.+DK1*DK1/(PI*PI-DK1*DK1) 

c*** 

TEST  FOR  A  STRING 

0150 

DATA  (I,J) =-FG/(SLAP*SLAP  +  DK1*DK1) 

0151 

GO  TO  1 

0152 

201 

CALL  CALCUL  (0) 

0153 

1 

CONTINUE 
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0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 


0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 

0181 

0182 

0183 

0184 

0185 

0186 

0187 
0188 
0189 
0190 
0191 
0192 
0193 
0  194 
0195 
0196 


IF  (NTEST.NE.1)  GO  TO  202 
WRITE  (6,507) 

CALL  FOURT  (DATA, NN, 2, 1, 1,0) 

DO  15  KI=  1,N2 
B=  -SI*PI/N*  (N-1)*  (KI-1) 

DO  15  KJ=  1,M 

C=  -SI  +  PI/H*  (M- 1)  *  (KJ-1) 

T=  PI/T  K*  ( K  J-1 ) 

DATA  (K1,KJ)  =  DATA  (KI,KJ)  *XK*TK/PI/PI/N/M*CEXP  (B)  *CEXP(C) 

15  DATA  (KI ,  KJ)  =  DATA  (KI, KJ) *EXP (CLAP*T) 

WRITE  (6,506)  (  (DATA  (I,J)  ,  J=  1,M),  1=  1,N2) 

GO  TO  203 

C* *************************  ** 

C***  MAIN  LOOP  FOR  DIFFERENT  DEPTHS  X3 
202  DO  6  K3=  1,NK3 
X3=  (K3-1)*DX3 

C***  LOOP  FOR  CALCULATING  DISPLACEMENTS  AND  STRESSES  AT  A  GIVEN  DEPTH 
DO  4  K=  1 ,NP 
WRITE  (6,680) 

WRITE  (6,501)  X  3 
REWIND  2 
DO  7  1=  1 , N 
DO  7  J=  1,H 

READ  (2)  P1,P2,C1,C2,S1,S2,DK1,DK2 
B=  CEXP (-X3*P  1) 

C=  CEXP  (-X3*P2) 

C***  DISPLACEMENTS 

IF  (K.EQ.  1)  DATA  (I ,  J)  =  Cl*B+C2*C 
IF  (K.EQ. 2)  DATA  (I,J) =  C  1*S  1*B+C2*S2*C 
C***  STRESSES 

IF  (K.EQ.  3) 

1  DATA  (I,J) =  (Cl 3*SI*DK 1-C33*P1*S1) *Cl*B> (Cl 3*SI*DK1-C3 3*P2*S2) *C2*C 
IF  (K.EQ. 4) 

1DATA  (I,J)=  (Cl  1  *S I*DK  1  -C13*P1*S1)  *C  1*B+  (C 1 1* SI*DK  1-C  1  3* P2*S 2)  *C2*C 
IF  (K.EQ. 5) 

1  DATA  (I,  J)  =  C55*  (  (~P1+SI*DK1*S1)  *Cl*B+(-P2  +  S_T<cDKl*S2)  *C2*C) 

IF  (K.I.T.6)  GO  TO  7 

C***  DISPL.  FOR  A  BEAM  ON  AN  ELASTIC  FOUNDATION 
EF=0 .  1 

C***  FORCING  FUNCTION 

FG=  FI*FG1* (1.+CEXP(-SI*DK2*T0) ) / (PI*PI-DK2*DK2*TC*T0 ) * 

14./  (DK 1*A0) **2*  (CSIN  (DK1*A0) /DK 1-AO*CCOS ( DK1* AO) ) 
D=W*CK1*DK1/12.0*(ES*W*W*DK1*DK1-RO*W*W*DK2*DK2)-RO*W*DK2*DK2 

1  +EF 

IF  (K.EQ.6) 

1DATA  (I,  J)  =FG/D 
7  CONTINUE 

CALL  FOURT (DATA ,N N, 2, 1 , 1 , 0) 

DO  3  KI=  1  ,N2 
B=  -SI*PI/N*  (N-1)  *  (KI-1) 

DO  3  KJ=  1 ,  M 
C=  -SI*PI/M* (M-1) * (KJ-1) 

T=  PI/TK*  (K J-  1 ) 

DATA  (KI , K J) =  DATA (KI,K J) *XK*TK/PI/PI/M/N*CEXP (B) *CEXP(C) 

E=  CL AP*T 

3  DAT  A  (KI,KJ)  =  DATA  (KI,  KJ)  *EX  P  ( E) 
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IV  G  LEVEL 

21  HAIN  DATE 

c*** 

CALCULATE  THE  PHOTOELASTIC  FRINGE  ORDER 

0197 

DO  290  1=1,40 

0198 

DO  290  J= 1, 32 

0199 

IF  (K.EQ.3)  T33 (1,3) =  REAL  (DATA (1,3) ) 

0200 

IF  (K.EC.4)  Til  (1,3)  =REAL  (DATA  (I,J)  ) 

0201 

IF  (K.EQ.5)  T13(I,J)=REAL(DATA(I,J)) 

0202 

IF  (K.EQ.5)  GO  TO  280 

0203 

GO  TO  290 

0204 

280 

CONTINUE 

0205 

FNGE=  (Til  (I,J)-T33(I,J))**2  +  4.0*T13(I,J)*T 

0206 

FRNGE  (I,J)  =SQBT  (FNGE) 

0207 

290 

CONTINUE 

0208 

214 

WRITE  (6,511)  K 

c*** 

FIND  THE  HAXIHUM  VALUE 

0209 

RS=  1.E-3 

0210 

DC  14  I  =  .1, 40 

0211 

DO  14  J=  1,32 

0212 

S=  DATA  (I,J) 

0213 

TP=  REAL  (S) /RS 

0214 

IF  (ABS  (TP)  .LT.  1  .)  GO  TO  14 

0215 

RS=  REAL  (S) 

0216 

14 

CONTINUE 

0217 

WRITE  (6,516)  RS 

C 

PRINT  REAL  PART  OF  DISPL.  AND 

0218 

NIJ=N/4 

0219 

NIJ=10 

0  220 

DC  310  1=1, NIJ 

0221 

11=1-1 

0222 

WRITE  (6,600)  11 

0223 

DO  300  J= 1 , 32 

0224 

SR=REAL  (DATA  (I,J)  )/BS 

0225 

SIdG=AIM AG (DATA (I, J) )/RS 

0226 

RDATA  (J)  =SR 

0227 

CIDATA(J)  =  SIHG 

0228 

300 

CONTINUE 

0229 

W  RITE  (6,620)  (BE  AT  A  (L)  ,L=  1 ,  M) 

0230 

WRITE  (6,650) 

0231 

WRITE  (6,620)  (CIDATA  (L)  ,  L=1 , 1 1 ) 

0232 

310 

CONTINUE 

0233 

WRITE  (6,680) 

0234 

IF  (K.EQ.1)  WRITE  (6,630) 

0235 

IF  (K.EQ.2)  WRITE  (6,631) 

0236 

IF  (K.EQ.3)  WRITE  (6,6  32) 

0237 

IF  (K.EQ.4)  WRITE (6,6  33) 

0238 

IF  (K.EQ.5)  WRITE (6,634) 

0239 

IF  (K.EQ.6)  WRITE  (6,635) 

0240 

WRITE (6,640)  DX, DT 

c*** 

PLOT  THE  RELATIVE  VALUES 

0241 

DO  12  J=  1,N 

0242 

DO  13  1=1,40 

C243 

S=  DATA  (I,  J) 

0244 

13 

MH(I)  =  REAL  (S)  /RS*  100 

0245 

WRITE  (6,515)  MM 

0246 

12 

CONTINUE 

c*** 

PLOT  A  HAP  OF  PHOTCELASTIC  FRINGE  ORDER 

0247 

IF  (K.LT.5)  GO  TO  4 
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0248 
0249 
0250 
0251 
0252 
0253 
0254 
0255 
0256 
0257 
0258 
0259 
0260 
0261 
0262 
0263 
0  264 
0265 
0266 
0267 
0268 
0269 
0270 

0271 

0272 
0273 
0274 
0275 
0276 
0277 
0278 
0279 
0  280 
0281 
0282 
0283 
0284 
0285 
0286 
0287 
0288 
0289 

0290 

0  291 
0292 
0293 


0294 

0295 

0296 

0297 


IV  G  LEVEL  21 


MAIN  DATE  3  74186  11/28/2 


RS= 1 . E-3 
DO  312  1=1,20 
DO  312  J=  1 , 20 
SS=FRNGE (I, J) 

TP=SS/FS 

IF (AES (TP) . LT. 1 . )  GO  TO  312 
RS=SS 

312  CONTINUE 

WRITE  (6,680) 

WRITE  (6,516)  RS 
WRITE (6,690) 

DO  320  J= 1 , 32 
DO  3  15  1=1,40 
SS=FRNGE (I,J) 

315  MM  (I) =SS/RS* 100 
WRITE  (6,515)  MM 

320  CONTINUE 
4  CONTINUE 
6  CONTINUE 

100  FORMAT  (8E10.4) 

101  FORMAT  (2I5,3F10. 3) 

102  FORMAT  (315) 

490  FORMAT (5X, 55H  ELASTIC  CONSTANTS  C 1 1 ,C22 ,C 33 ,C44, C55 , C66 , C 12 , C 1 3 , C3 
12 ,/9E12. 4 ,/) 

491  FORMAT (5X, 25H  DENSITY  OF  PLATE  GH/CC  — , FI  0. 4, 5X, 20H  FIBER  LAYUP  AN 
1GLE — , FI  0 . 4) 

500  FORMAT  (2F 1 2. 4, 4E1 5. 4) 

501  FORMAT  (5H  X3=  F10.4) 

502  FORMAT  ( 1 0 H  STRESSES  I3/(8E15.7)) 

503  FORMAT  (5X5H  DATA,  15X10 HCONSTANTS  , 20X 10HPARAMETEBS  /8E15.7) 

504  FORMAT  (5X8HNSTRESS=  I3,5X4RNX3=  I4,5X4HDX3=  F10.4,/) 

505  FORMAT  (14H  SIMPLE  POLES  ) 

506  FORMAT  (8E15.7) 

507  FORMAT  (15 H  THIS  IS  A  TEST  ) 

5C8  FORMAT  (30H  LAPLACE  INVERSION  PARAMETER3  F12.4) 

509  FORMAT  (28H  SECCNE  ORDER  BRANCH  POINTS  ) 

510  FORMAT (5X, 27H  LONG. WAVE  SPEED  IN  BEAM  =  ,F12.3,7H  CM/SEC) 

511  FORMAT  ( 14H  DISPLACEMENTS  14) 

512  FORMAT  (24H  RAYLEIGH  SPEED  CS/CR*I  ) 

513  FORMAT  (8F15.8) 

514  FORMAT  (/50X  10HWITH  STRIP) 

515  FORMAT  (2X, 4013) 

516  FORMAT  (16H  MAXIMUM  VALUE3  E15.7) 

520  FORMAT (/, 5X,23H  IMPACT  PRESSURE (PSI) =  , El 2 . 4 , 5X , 25H  1/2  CONTACT  LE 
1 NGTH  (CM) =  ,F12. 4,5X,23H  IMPACT  DURATION (SEC) =  ,E12.4,/) 

521  FORMAT (/,5X, 18H  DENSITY  OF  BEAM=  ,F12.4,5X, 17H  NORM. THICKNESS3  , 
1F12.4,5X,15H  BEAM  MODULUS3  ,E12.4) 

600  FORMAT (/,  20X, 1 8HNORM ALIZED  DIST.3  ,14) 

620  FORMAT  (11F10.5) 

621  FORMAT (5X, 20H  W A VELEN (TIME-S EC) 3  ,E12.4 ,5X,20H  WA VELE N (DI ST-CM  )  = 
1,F12.4,/,5X,23H  MAX  FREQ  NO(NON  DIM) 3  , Fl 2. 4 , 5X, 2 3H  MAX  WAVE  NO  (NC 
2N  DIM)3  ,F12.4,/) 

630  FORMAT  (30 X,16H  DISPLACEMENT  Ul,/) 

631  FORMAT  (3CX, 16H  DISPLACEMENT  U3,/) 

632  FORMAT  (30X,1 1H  STRESS  T33,/) 

633  FORMAT  (30X,1 1H  STRESS  Til,/) 


42 


OUTRAN  IV  G  LEVEL  21 


MAIN 


DATE  -  74186 


11/28/2 


C298 

0299 

0300 

0  30  1 
0302 
0  30  3 

0304 

0305 

0306 


634  FORMAT  (30X,1 1H  STRESS  T13,/) 

635  FORMAT (30X,41H  TEST-DISPL . OF  BEAM  ON  ELASTIC  FODNDATION ,/) 

640  FORMAT (20X, 14H - »  XI,  DX=  ,F1 2 . 4, /, 20X , 2H  |,/,20X,2H  |,/,2QX, 

1  2H  |,/,20X,2H  V,/,20X,2H  V,/,20X,9HTIME,DT=  ,E12.4,//) 

650  FORMAT ( 10H  IMAG  PART) 

680  FORMAT  (1 H 1 )  „  ,,  „  ,  _ 

690  FORMAT  (/, 20X,  17H  FRINGE  ORDER  MAP,/ *  20X , 3  1H  SQRT  (  (T1 1  -T33)  **2+4  .  -j* 

1T13*T13) ) 

710  FORMAT (  5X,6H  VEL=  ,F12.4,5X,5H  T0=  ,F12.4,5X,4H  A=  ,F12.4,5X, 

15 H  F0=  ,F12. 4,/) 

203  STOP 
END 
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COC-1 
0002 
000  3 
0004 
0005 
0006 
0007 
0008 
0009 
00  10 
0011 
0012 

00  13 
0014 
0015 
0016 
0017 
0018 
0019 
0020 
0021 
0022 
0023 
0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 
0034 
0035 
0036 
0037 
00  38 
00  3  9 
00  40 
004  1 
0042 
0043 
0044 
0045 
C046 
0047 
004  8 

0049 


IV  G  LEVEL  21  MAIN  DATE  ~  74186  11/28/2 


SUBBOOTINE  CALCUL  (NPOLE) 

COMMON  /MMC/  D, DK1 

COMMON  /MC/DK2,FGl,T0,Cl1 ,C 1 3,C33 ,C 55 , RHO ,R0, H, ES , NSTRI P , AO 
COMPLEX  B 1 , B2 

COMPLEX  G0,GG1, HO, HOI , G 10 ,G  1 1 ,H1 0,H 1 1 
COMPLEX  G,H,G1,H1 

COMPLEX  DK1,DK2,P1,P2,S1,S2,C1,C2,SI,B,C,D  ,FG 
COMPLEX* 16  BS 
PI-  3.14159265 
SI  =  CMPLX  (0 . ,  1 . ) 

A=  C33*C55 

B=  C33* (RHO*DK2*DK2-C11*DK1*DK1) +C55* (P HO*DK2*DK2-C55 *DK 1 *DK 1)  + 
1DK1*DR1*  (C13+C5  5) **2 

C=  (RH0*DK2*DK2-C11*DK1*DR1) * (RHO*DK2 *DK2-C55*DK1 *DK 1 ) 

BS=  1.D0*B*B-1 .D0*4.*A*C 
D=  CESQRT(BS) 

Pl=  . 5/A*  (-B*D) 

P2=  C/PI/A 
P1=  CSQRT(PI) 

IF  (REAL(PI)  . LT. 0.0. AND. NPOLE. NE.  2)  Pl  =  -Pi 
P2=  CSQRT  (P2) 

IF  (REAL (P2) .LT. 0.0. AND. NPOLE. NE. 2)  P2=-P2 

S1  =  -SI/DK1  /Pi  /  (C13+C55) *  (RH0*DR2*DK2-C11*DK1*DR1+C55*P1*P1)  ‘ 

S2=  -SI/DK1  /P2/  (C13+C55) * ( RHO*DK2*DK 2-C 1 1*DK 1*DK1 +C55*P 2*P2) 

IF  (NPOLE. EQ. 2)  GO  TO  201 
C***  FORCING  FUNCTION 

FG=  PI*FG1 *  (1. +CEXP (-SI*DK2*T0) ) / (PI*PI-DK2*DK2*T0*T0) * 

14./  (DK 1  *A0 )  **2*  (CSIN  (DK1  *A  0)  /D.K  1  -  A0*CCOS  (DK  1  *  AO)  ) 

201  IF  (NSTRIP.EQ. 1)  GO  TO  202 

D=  (SI*CK  1*S1-P 1) *  (SI*DK1*C13-P2*C33*S2) - (SI*DK1*S2-P2)  * (SI*DK1* 
1C13-P1*C33*S1) 

IF  (NPOLE. EQ. 2)  RETURN 
206  Cl=  FG/D* (5I*DK1*S2-P2) 

C2=  FG/D* (P1-SI*DK1*S1) 

GO  TC  203 
202  CONTINUE 

C***  MATRIX  FOR  STRIP  PROBLEM 

B1=W*H*(ES*DK1*DK1-R0*DK2*DK2) *DK 1*DK 1/1 2 . 0- RO*DK2*DK 2 
E 1= H*B  1 

B2=-N* (ES*DK1*DK1-R0*DK2*DK2) 

G 10=  C55*(SI*DK1*S1-P1) 

H 1 0=  C55*  (SI*DK1*S2-P2) 

G11=E2*(1.0+SI*DKl*S1*W/2) 

H 1 1 =E2*  (1.0*SI*DKl*S2*N/2) 

G0=- (SI*DK1*C13-P1*S1*C33) 

H0=- (SI*DK1*C13-P2*S2*C33) 

G01=-SI*DK1*M*G1 0/2 . 0  +  Bl*Sl 

HO 1=-SI*DK1*N*H 10/2.0  +  B1*S2 

H 1  =  H  10+H 1 1 

G1=G10*G11 

G=GO+G01 

H=H0+H0 1 

D=G1C*H0-G0*H10  +  (H0*G11  +  H01*G10-G0*H1 1-G01*H10)  + ( H 0  T* G 1 1 - H 1 1*G01 

D 

IF  (NPOLE. EQ. 2)  RETURN 


FORTRAN 

0050 

0051 

0052 

0053 

0054 
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IV  6  LEVEL  21  CALCOL  DATE  =  74186 

C1=  ~FG/D*H 1 
C2=  FG/D*G1 

203  WRITE  (2)  Pi ,P2 ,C 1 ,C2 ,S 1 , S2, DK1, DK2 
BETUFN 
END 
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FIGURE  2b 

Comparison  of  In-plane  plate  edge  wave  speed  with  body  wave  speeds  versus 
fiber  layup  angle  for  55%  graphite  fiber/ epoxy  matrix  composite. 
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FIGURE  4 

Edge  Stress  t„  Surface  Wave  for  0°  Fiber  layup  Angle:  55% 
graphite  fiber/ epoxy  matrix  composite 


FIGURE  5 

Edge  Stress  Surface  Wave  for  +_  15°  Fiber  Layup  Angle:  55%  Graphite  fiber/expoxy 

matrix  composite. 
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FIGURE  6 


Computer  Map  of  Edge  Impact  Pressure  t  Impact  Time  87  ysec. 
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Computer  Map  of  Edge  Stress  tn  Surface  Wave  in  the  Space-Time  (x, ,  t)  Domain 
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FIGURE  9 

Stress  t  versus  Time  at  Different  Distances  from  the  Edge 
under  the  Contact  Point. 


FIGURE  10 

Maximum  Stresses  versus  Layup  angle  for  55%  Graphite  Fiber/epoxy 
Matrix  Composite. 
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FIGURE  11 


Distribution  of  Stress  t  Along  the  Edge  for  Fiber  Layup  angle* 
0°,  +  30,  +  45°.  ** 


NORMALIZED  STRESSES  IN  COMPOSITE 
AT  STRIP -EDGE  INTERFACE 


ALIZED  STRESSES  IN  COMPOSITE 
r  STRIP -EDGE  INTERFACE 


DATA  FDR  55%  GRAPHITE  FIBER/ EPOXY  MATRIX 
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FIGURE  12b 

Effect  of  Edge  Strip  Thickness  on  Interface  Stresses. 


0°  LAYUP  ANGLE  ALUMINUM  STRIP 
X3=  0  ,  a  =  I  cm  ,  T0  =  35 /a sec 


FIGURE  13 

Effect  of  Edge  Strip  Thickness  on  Interface  Stresses. 


FIGURE  15 

Effect  of  Edge  Strip  Thickness  on  the  Rayleigh  Edge  Wave  Shape. 


